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In  the  framework  of  the  Comprehensive  Test  Ban  Treaty,  development  of  reliable 
methods  to  discriminate  between  underground  nuclear  explosions  and  earthquakes 
at  regional  distances  (less  than  2500km)  continues  to  be  very  important  especially 
in  connection  with  the  last  (in  May,  1998)  nuclear  explosions  conducted  at  Indian 
and  Pakistan  test  sites.  Since  the  lithosphere  is  a fractal,  we  suppose  the  signals, 
which  propagate  through  the  media,  inherit  its  ’self-similar’  (scaling)  features.  We 
assumed  that  these  features  of  explosions  and  earthquakes  or  their  topological 
reconstructions  (embeddings)  have  to  be  different.  Scaling  reflects  correlations  of 
more  high  order  then  it  is  possible  to  estimate  by  linear  discriminating  methods  and 
can  be  used  as  base  of  non-linear  discrimination.  We  propose  to  build  a universal 
geometrical  model  of  a seismic  signal  using  the  canon  algorithm  of  F.  Takens  and 
to  estimate  scaling  of  the  model.  The  scaling  features  were  used  as  patterns  of 
seismic  signals  for  entering  them  into  an  artificial  neural  network.  Records  of 
nuclear  explosions  and  earthquakes  from  different  regions  were  included  into  the 
training  set.  The  net  was  trained  to  classify  types  of  seismic  events.  Results 
have  shown  89%  correct  classification  of  the  unknown  signals.  As  additional  tools 
for  distinguishing  between  nuclear  explosions  and  earthquakes  we  propose  to  use 
Hurst’s  method  and  the  cross  correlation  method.  Results  of  using  these  methods 
are  demonstrated  on  examples  of  some  explosions  and  earthquakes. 


1 Introduction 


The  nuclear  explosion  discrimination  problem  continues  to  be  very  important  espe- 
cially in  connection  with  the  last  ( in  May,  1998)  nuclear  explosions  conducted  at 
Indian  and  Pakistan  test  sites.  Existing  regional  methods  of  seismic  event  discrim- 
ination 1’2<3'4’5  are  based  on  the  comparative  analysis  of  spectral  characteristics  of 
two  main  components  (P-wave  and  S-wave)  of  a seismogram  (see  Fig.  1).  How- 
ever, these  parameters  are  very  sensitive  to  non- uniformity  of  the  lithosphere  and 
the  astenosphere  and  depend  on  the  location  of  the  event  and  the  path  of  a signal 
propagation.  Moreover,  modern  technology  of  nuclear  testing  complicates  distin- 
guishing between  nuclear  explosions  and  earthquakes. 

It  is  known  that  the  lithosphere  exibites  fractal  features  6,7  in  a wide  range  of 
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Figure  1.  A seismogram  of  the  Chinese  earthquake,  24.07.86 


scales  - from  parts  of  millimetre  of  geological  grain  to  tens  thousands  kilometres 
of  a tectonic  plateau  and  has  hierarchy  of  self-similar  blocks.  Response  of  the 
fractal  lithosphere  depends  on  source’s  geometry  (radiation  patterns)  of  a seismic 
event  more  than  on  its  yields.  Moreover,  an  earthquake  centre  is  formed  in  the 
media  during  the  long  time.  Features  of  the  media  are  changed  by  the  moment  of 
the  earthquake,  which  occurs  in  the  preliminary  stressed  media.  An  explosion  is  an 
artifact  for  the  media  and  is  not  connected  with  any  preliminary  changes  of  external 
parameters  of  the  media.  Since  dynamic  scenarios  of  explosions  and  earthquakes 
are  different,  we  assumed  that  scaling  features  of  seismic  signals  or  their  topological 
reconstructions  (embeddings)  have  to  be  different.  Scaling  reflects  correlations  of 
higher  order  than  it  is  possible  to  estimate  by  linear  discriminating  methods,  and 
can  be  used  as  base  of  non-linear  discrimination.  We  propose  to  build  a universal 
geometrical  model  of  a seismic  signal  using  the  canon  Takens’  algorithm  8’9  and  to 
use  its  scaling  features  as  attributes  of  seismic  event  discrimination. 

The  structure  of  this  paper  is  as  follows.  In  Sec.  2 we  study  ID  scaling  features 
of  seismograms.  In  Sec.  3 we  study  correlation  dimension  of  reconstructions  of 
explosions  and  earthquakes.  In  Sec.  4 we  use  the  scaling  features  of  embeddings 
of  seismic  signals  as  patterns  for  a neural  network.  In  Sec.  5 we  describe  the 
cross  correlation  method  as  a tool  for  seismic  discrimination.  Our  findings  are 
summarized  in  the  conclusion  section. 
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2 Hurst’s  exponents  of  seismograms 

We  started  from  the  study  of  self-similar  properties  of  the  processes  in  seismic 
sources  applying  Hurst’s  method  10  to  the  seismograms  of  nuclear  explosions  and 
earthquakes.  Let  us  remind  the  classical  notion  of  self-similarity  for  a random 
process11.  X(t ) is  a self-similar  process  in  Sft  if  there  exists  a sequence  of  positive 
real  numbers  cr  such  that,  for  any  r > 0, 


X(t)  = CrX(rt),t  £ 1ft  (1) 

where  = denotes  equality  of  all  finite-dimensional  distribution.  This  equation  is 
a statement  of  invariance  of  X(t)  under  the  group  of  affine  transformations  X —> 
crX,t  —>  rt,Cr  > 0.  Since  ca{,  = cact  and  cj  - 1,  cr  must  have  the  form  rH  for 
some  H > 0 and  this  formula  might  be  modified  as 

X(t)  = r~HX(rt),t£  SR  (2) 

where  H is  the  Hurst  exponent.  Traditionally  it  is  estimated  by  the  rescaled  range 
method  10 . 

We  analyzed  seismograms  of  the  nuclear  explosions  conducted  in  India  and  Pak- 
istan and  Tibetian  earthquakes.  Fig.  2 represents  the  Hurst’s  exponents  for  different 
types  of  seismic  events.  Analysis  of  the  Hurst’s  exponents  allowed  to  conclude  that 
graphs  of  the  earthquakes  are  more  regular  than  graphs  of  the  explosions  because 
H estimated  for  earthquakes  equals  1 on  larger  range  of  scales  (i.e.  logT  £ [0, 1.5]) 
than  for  explosions  (i.e.  logT  € [0, 1]). 

It  is  necessary  to  note  that  the  procedure  described  above  is  suitable  only  for 
records  with  low  level  of  noise. 

3 Correlation  integral  as  a tool  for  distinguishing  between  different 
seismic  sources 

The  estimation  of  the  correlation  dimension  12  of  attractors  reconstructed  directly 
from  experimental  time  series  is  often  used  means  of  gaining  information  about 
the  nature  of  the  underlying  dynamics.  It  is  proposed  that  a scalar  time  series 
y(t)  = = 1,2,  ...,1V  have  been  generated  by  a smooth  dynamical  system 

x(f  + t)  = fr(x(f))  defined  on  a manifold  M with  an  attractor  A £ M such  that 
yi  = /i(x(it))  where  h : M — » 1ft  is  Lipschiz  function.  Then  the  reconstruction 

y(i)  {Uii  Vi+ 1>  •••>  Vi+m—  l}  = 

= {Mx(*))>  Mx(*  - t)),  h(x(t  - 2 t)),  ...,  /i(x(t  - (m  - 1)t))} 

defines  a delay- coordinate  map  A : M — > and  A — > is  the  reconstructed 

attractor.  Takens’  theorem  8,9  ensures  that  if  m > 2D  where  D is  the  box-counting 
dimension  (or  capacity)  of  A,  then  the  map  A is  embedding,  i.e.  is  one-to-one,  and 
also  an  immersion  with  a precision  upto  assumption  about  prevalence  13 . 

For  Takens’  theorem  to  be  valid,  we  need  to  assume  that  both  the  dynamics  and 
the  observations  are  autonomous  (so  that  f and  h depend  on  x only).  Unfortunately, 
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Figure  2.  Hurst  exponent  H curves  for  different  seismic  events:  Pakl  - nuclear  explosion  30.05.98 
(Pakistan),  ind  - nuclear  explosion  11.05.98  (India),  tibl  - earthquake  (Tibet),  tib2  - earthquake 
(Tibet),  tib3  - earthquake  (Tibet) 


these  assumptions  fail  to  hold  for  seismic  sources.  However,  there  exists  Takens 
embedding  theorem  for  forced  and  stochastic  systems  14 . We  applied  the  technique 
of  reconstruction  described  above  basing  on  the  last  theorem  and  assuming  that 
the  seismograms  are  typical  and  contain  all  the  character  scales  of  the  dynamics  of 
the  seismic  sources.  It  means  that  in  some  sense  the  attractor  of  the  seismic  source 
exists. 

We  reconstructed  some  old  non-camouflaged  nuclear  explosions  and  found  out 
that  their  embeddings  in  3?2  visually  differ  from  ones  of  earthquakes  (see  Figures  3- 
4).  Unfortunately,  this  difference  for  records  of  last  explosions  (Indian  and  Pakistani 
explosions)  is  not  visible  (see,  for  example,  Fig.  5).  Therefore,  it  is  reasonable  to 
use  numerous  characteristics  of  these  reconstructions.  The  most  popular  tool  for 
description  of  the  embedding  obtained  is  the  correlation  integral  12 . 

Let  y (i)  £ 9?m  be  a point  on  the  attractor  Ar.  The  correlation  integral  is 
defined  as  proportion  of  pairs  of  points  of  no  more  than  distance  e apart.  That  is, 

Cm(e)  = ^ 0(lly(i)  - y(fc)ll  - g)  (4) 

' ' 3 k 

Here:  © - Heaviside  step  function,  the  symbol  ||  * ||  always  denotes  ’sup’  norm 
on  3ft771. 

There  is  scaling: 


Cm(e)  oc  ev 


(5) 


Figure  3.  Embedding  in  5R2,  nuclear  explosion,  Figure  4.  Embedding  in  3t  , earthquake, 
STS,  25.04.82.  China,  24.06.86. 


Figure  5.  Embedding  in  5R2,  the  Indian  nuclear  explosion,  recorded  by  NIL  station,  30.05.98. 


for  e — > 0.  The  slope  of  the  correlation  integral  is  called  correlation  dimension 


A log  C(e) 

v = lim  — — (6) 

e->0  Alog£ 

The  estimated  v typically  increases  with  m and  reaches  a plateau  (in  the  best  case) , 
on  which  the  v estimate  becomes  relative  constant. 

It  is  expected  that  because  of  the  spherical  symmetry  and  small  sizes  of  sources 
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of  underground  explosions  the  main  part  of  seismic  energy  is  contained  in  compres- 
sional  P-waves  having  high  frequency  content  in  comparison  with  earthquakes.  At 
the  same  time  earthquakes  differ  by  generating  intensive  shear  S-waves.  These  dif- 
ferences are  traditionally  used  for  choosing  diagnostic  parameters  for  seismic  event 
discrimination.  Hence,  we  studied  the  signals  for  both  types  of  waves  separately. 

We  processed  records  of  underground  nuclear  explosions  conducted  at  Semi- 
palatinsk  Test  Site  (STS, Kazakhstan)  and  Lop  Nor  (China)  and  earthquakes  in 
China  and  Altay  (Russia).  The  seismograms  were  recorded  by  Kazakhstani  seis- 
mic stations  BRVK,  VOS,  ZRN,  CHK  located  in  the  North-West  and  TLG  - in 
the  South  of  Kazakhstan.  In  addition,  seimograms  of  NIL  station  (Pakistan)  were 
processed. 

Our  experience  showed  that  correlation  integrals  for  P-wave  reconstructions 
are  more  informative  than  ones  for  S-waves:  analyzing  correlation  integrals  for  S- 
waves  we  didn’t  find  out  notable  features  discriminating  different  seismic  sources. 
Fig.  7 shows  the  typical  correlation  integrals  for  P-waves  of  nuclear  explosions  and 
earthquakes.  It  is  seen  that  the  correlation  dimensions  calculated  for  different  types 
of  seismic  events  are  different.  As  a rule,  the  correlation  dimensions  are  higher  for 
explosions  than  for  earthquakes. 


Figure  6.  Slopes  of  correlation  integrals  (m  = 
8 — 16,  r = 2),  P-wave  of  the  nuclear  explosion, 
Semipalatinsk  Test  Site,  25.04.82. 


Figure  7.  Slopes  of  correlation  integrals  (rn.  = 
7 — 14,  r = 2),  P-wave  of  the  Chinese  earth- 
quake, 24.07.86. 


In  general  case  the  slopes  exibit  complex  behavior  of  u versus  e Fig.  8.  For 
complex  systems  it  is  possible  that  more  than  one  scaling  exists  for  different  e. 
The  complexity  of  seismic  record  slopes  might  be  explained  by  presence  of  different 
kinds  folding  effects  in  3?m  17 . For  the  reasons  outlined  above  we  used  the  slope’s 
form  instead  of  the  value  of  the  correlation  dimension,  i.e.  the  function  u = r'(loge). 

We  used  this  function,  as  discriminating  attribute  of  seismic  signals  for  training 
an  artificial  neural  network  to  classify  seismic  sources. 
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Figure  8.  Slopes  of  correlation  integrals  (m  = 10  — 25,  r = 2),  P-wave  of  the  Indian  nuclear 
explosion  recorded  by  NIL  station  30.05.98. 

4 Neural  network  using 

In  the  last  decade  artificial  neural  networks  (ANN)  have  become  the  most  popular 
tool  for  pattern  recognition  tasks  19,20  in  general,  and  in  seismic  sources  distin- 
guishing in  particular  3,4,s.  One  can  consider  a neural  network  as  a distributed 
dynamical  system  consisting  of  a set  of  non-linear  processing  elements  (formal  neu- 
rons) connected  according  to  a certain  architecture.  The  formal  neuron  is  able  to 
pick  up  input  vectors  x = {xj},  j = 1,  2, ...,  N,  estimate  their  scalar  multiplications 
with  weights  w = {m*},  i = 1,2, M: 

S = YlxiWi  (7) 

and  transform  S in  accordance  with  an  activation  function  </>  into  the  output  vec- 
tor of  the  neuron  y — cj)(S).  Such  a neural  network  might  be  trained  to  recognize 
unknown  patterns.  Network  training  process  is  based  on  adjusting  the  weight  con- 
nections between  related  values  of  inputs  and  targets  (desirable  outputs)  so  as  to 
minimize  an  error  function  19 . A set  of  input-targets  values  is  called  a training  one. 
There  exist  two  main  types  of  ANN:  fully-connected  nets  and  perceptrons  20 . 

In  order  to  check  the  algorithm  described  in  the  previous  section  we  trained 
the  fully-connected  artificial  neural  network  ’’MultiNeuron”  developed  by  Russian 
scientists  20  using  60  examples  of  the  correlation  integral  slopes  of  seismic  records 
mentioned  above  as  input  patterns.  The  net  was  tested  on  10  seismograms  of 
different  underground  nuclear  explosions  and  earthquakes  which  were  not  included 
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into  the  training  set.  Results  of  testing  showed  89%  of  the  correct  classification  of 
the  signals. 

5 Cross-correlation  sums  as  a tool  for  distinguishing  between 

different  seismic  sources 

In  this  section  we  represent  results  of  our  study  of  non-linear  cross  correlation  sums, 
which  are  the  generalized  correlation  integrals  and  are  used  here  for  estimating 
similarity  of  two  probabilistic  measures. 

Let  {ui),  {zi},  i = 1,2, ...,  N denote  two  different  observed  time  series.  Internal 
dynamics  of  systems  in  a phase  space  M which  produces  these  series  as  typical 
mapping  M — » 5J  is  unknown.  However,  it  is  proposed  that  there  exist  chaotic  or 
quasiperiodic  attractors  of  these  systems,  dimensions  of  which  dy  and  dz  are  low 
enough  for  applying  the  reconstruction’s  methods. 

Let  and  A be  their  embeddings  into  3in,n  > 2 (dy  + dz).  Let  y and  z be 
corresponding  delay-vectors  selected  randomly  according  to  two  different  measures 
p and  p.  The  cross  correlation  integrals  are  defined  by15  : 

Cmp(£)  = J My)  J dp{z)&(£  - ||y  - z||)  (8) 

The  cross  correlation  sum  is  defined  similarly  by 

Cyz(£)  = AT-2^^0(£  - ||y-z||)  (9) 

y z 

Almost  surely,  Cyz  (£)  ~ > C>p(e)  for  N — > oo,  and  Cyz(e)  is  an  unbiased  estima- 
tor of  For  sufficiently  small  e and  for  absolutely  continues  measures  p and 

p one  can  show  15 

Clp(e)<C^(s)Cpp(s)  (10) 

In  practice,  this  inequality  is  used  to  calculate  the  cross  correlation  ratio  as  a 
similarity  measure 


r(e)  = 


Cyz(£) 

JCyy(E)CZZ(E) 


We  used  it  in  the  form  16 


(11) 


TUfii  ~ z(i)ll2Q(£  - lly(»)  - y(j)ll) 

Ei^0(£-  IlyOO  -y(j)ll) 

If  {yi}  and  {z2-}  are  related  by  identical  dynamical  scenarios  than  one  can  expect 
that  ||y(i)  — y(j)||  < £ =>  ||z (i)  - z(j)\\  ss  e.  If  no,  ||z(i)  -z(j)||  « Az,  where  Az  is 
the  average  distance  between  points  of  the  attractor  A ^ and  A'”’  does  not  depend 
on  e. 

Fig.  9 represents  the  cross  corelations  of  various  seismic  events.  The  cross  corre- 
lation between  the  Tibetian  earthquake  (09.01.90)  and  the  Indian  nuclear  explosion 
(11.05.98)  is  absent:  FT™  (e)  curve  is  practically  in  parallel  with  the  horizontal  axis. 
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Figure  9.  Cross  correlation  between  recon- 
structions of  the  Tibetan  earthquake  and  the 
Indian  nuclear  explosion. 


Figure  10.  Cross  correlation  between  recon- 
structions of  the  Indian  and  the  Pakistani  nu- 
clear explosions. 


The  Figure  10  demonstrates  strong  cross  correlation  between  the  Indian  explo- 
sion (11.05.98)  and  the  Pakistani  nuclear  explosion  (28.05.98).  In  this  case  we  can 
conclude  that  dynamical  scenarios  of  two  systems  observed  are  identical  ones. 


6 Conclusion 

We  have  introduced  non-linear  criteria  for  distinguishing  between  underground  nu- 
clear explosions  and  earthquakes  recorded  at  regional  distances  using  the  fractal 
approach. 

We  noted  that  behavior  of  Hurt’s  exponents  varies  for  graphs  of  different  seis- 
mograms in  dependence  on  scale:  the  regularity  range  for  the  earthquakes  is  larger 
than  for  the  explosions. 

We  found  out  that  embeddings  of  non-camouflaged  explosions  into  5ft2  visually 
differ  from  ones  of  earthquakes.  The  slopes  of  correlation  integrals  demonstrate 
complexity  of  scaling.  The  notable  plateau  is  absent  for  both  types  of  seismic 
sources  but  the  curves  of  the  slopes  are  different  for  explosions  and  earthquakes. 

The  correlation  dimension  curves  of  60  seismograms  of  different  seismic  events 
were  entered  into  the  neuroimitator  ’’MultiNeuron” , which  was  used  as  a classi- 
ficator  of  underground  nuclear  explosions  and  earthquakes.  The  testing  results 
obtained  have  shown  89%  correct  classification  of  10  seismic  events,  which  were  not 
included  into  the  training  set. 

Sometimes,  it  is  possible  to  identify  the  seismic  source  estimating  its  cross  cor- 
relation with  a test  example.  Analysis  of  relationships  of  attractors  in  the  phase 
space  allows  obtaining  additional  information  for  the  discrimination  task. 

Therefore,  our  experience  has  shown  that  the  non-linear  methods  could  be  suc- 
cessfully applied  to  regional  seismic  event  discrimination. 
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